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Abstract. We consider the effects of long-range temporal correlations in many- 
particle systems, focusing particularly on fluctuations about the typical behaviour. For 
a specific class of memory dependence we discuss the modification of the large deviation 
principle describing the probability of rare currents and show how superdiffusive 
behaviour can emerge. We illustrate the general framework with detailed calculations 
for a memory-dependent version of the totally asymmetric simple exclusion process as 
well as indicating connections to other recent work. 


1. Introduction 

Interacting particle systems in driven steady states are typically characterized by non¬ 
zero currents; among the recent advances in non-equilibrium statistical mechanics has 
been a considerable body of work on understanding the fluctuations of such currents. 
Indeed it is now well-established that, for Markovian dynamics, the probability of seeing 
a time-averaged current away from the mean is generically captured by a large deviation 
principle with “speed” t PEI [3]. However, models with some form of non-Markovian 
dynamics arguably describe better the long-range temporal correlations in many real 
scenarios PEIE]. In this direction, there is topical interest in both the typical behaviour 
and fluctuations for particle systems with memory. In particular, statistical physicists 
have recently studied a variety of memory-dependent random walkers in classical and 
quantum contexts, see e.g., PIHIEIIIQ] - some of these can be related to the reinforced 
random walks and Polya urn models found in earlier mathematical literature and 
reviewed, for instance, in [11]. Much less is known about non-Markovian many-particle 
systems but some aspects of the stationary-state properties (e.g., mean current as a 
function of density, conditions for a condensation transition) have been investigated for 
models with internal states or non-exponential waiting times [m ESI EH. Going beyond 
the typical behaviour, the current fluctuations in a temporally-correlated zero-range 
process have also recently been explored (and compared to the equivalent memoryless 
model) although exact analytical calculations proved possible only for a single site |15j . 
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In the present contribution we build on earlier work in [12] to show how an 
expansion about hxed points of the dynamics can yield valuable information about 
the fluctuations in a particular class of non-Markovian interacting particle systems, 
even when full solution appears a formidable task. Specihcally, this enables us to 
predict the speed of the current large deviation principle and hence the long-time scaling 
behaviour of fluctuations. We demonstrate this approach with perhaps one of the most 
famous models in non-equilibrium statistical mechanics: the totally asymmetric simple 
exclusion process (TASEP). Here we show how a current-dependent input rate leads to 
a modihed phase diagram including a super diffusive regime and we check our theoretical 
approximations against simulations and exact numerics. 

The remainder of the paper is structured as follows. In section [2] we introduce 
the framework of systems with current-dependent rates and indicate the connections 
to other recent works, including some of those mentioned above. In section [3] we 
perform a stability analysis of hxed points and make a Gaussian expansion to study 
the huctuations. The power of this approach is then illustrated by treatment of the 
TASEP in section H] before a concluding discussion and wider perspective in section |5l 
Finally, a short appendix provides a pedagogical treatment of a single-particle problem 
in order to demonstrate the formalism. 


2. Interacting particle systems with current-dependent rates 


We work within a discrete-space and continuous-time framework with the particle 
conhguration at time t labelled by a{t) and transition rate from state a to a' given 
by the matrix element Classical lattice-based many-particle models described in 

this way include exclusion processes (to which we will return later) [T7l fT8] . zero-range 
models [iHiiini, and inclusion processes m- 

In systems of this type, a time-integrated particle current J'{t) can be dehned as 
the net number of jumps across a given bond (or subset of bonds) from time zero up 
to time t. We choose the script style to indicate that the current is a functional of 
the stochastic history {(j(r),0 < r < t} but will suppress the explicit dependence on 
t where no confusion should arise. It is well-known that J' generically obeys a large 
deviation principle which can be loosely stated as 


Prob 





( 1 ) 


where ~ denotes logarithmic equivalence in the long-time limit. Iwij) is known as 
the rate function and t, somewhat misleadingly, referred to as the speed. The w 
subscript emphasizes that the rate function depends in some non-trivial way on the 
set of transition rates. Much recent industry has been devoted to calculating /^(j) 
for various models both within a “microscopic” lattice-based approach [2] and in the 
“macroscopic” hydrodynamic limit [3]. 

Here, following [16], we introduce an element of memory by considering a class 
of models in which the rates at time t depend on the current up to time t. To be 
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precise, the rates tCo-'.o- now depend on the time average J'/t of some specihed particle 
current and will be denoted by Obviously, the functional dependence on the 

current must be chosen so that the rates always remain positive. To avoid singularities 
at time zero we also assume that the time-averaged current starts at some fixed value 
jo at a time to which is small compared to the overall measurement time t. Physically, 
taking the condition 0 to ^ ^ excludes any initial transient behaviour which could be 
governed by different rates. A conceptually simple generalization is to the case where 
rates depend on multiple currents (e.g., currents measured separately across different 
bonds or in different directions) - much of the following analysis can be extended to that 
situation although we shall not consider it in detail. We emphasize that we specialize 
here to models with a functional dependence on the single variable Jjt {time-averaged 
current), rather than a more general dependence on jT" and t separately. The closely 
related scenario of feedback depending on the time-integrated current has also recently 
been explored, for example, in a quantum context [21]. 

For illustrative purposes we now consider a single particle, i.e., a random walker, 
with this type of memory and endeavour to describe its connection to a range of models 
in the literature which may, at first sight, appear rather disparate. The natural “current” 
for a particle on a one-dimensional lattice is just the net number of steps made in one 
direction, say towards the right, and we now assume left and right hopping rates which 
depend on the time average of this quantity (in other words, on the particle velocity). 
In the discrete-time version of this picture, the dependence is thus on the particle’s 
position divided by the number of time steps elapsed. Dynamics in this category 
includes the “elephant” random walker of [7] as well as several other recent random 
walk scenarios [22l |23l [24], the voting model of [25], and some aspects of the behaviour 
in a discrete-choice model with dependence on the peak of past experience [26]. In 
fact, mathematically, these are all essentially equivalent to the much older Polya urn 
problem izg in which the probability for picking a black or white ball depends on the 
relative number (fraction) of such balls chosen in the past. If the functional form of the 
dependence is non-linear, then one has a generalized Polya process, for overviews see, 
e.g., [II1|28]. In passing, we remark that such models can also be considered as a limiting 
case of binary Markov chains with memory of a finite number of steps; see e.g., [2^ [30] 
and note that the latter reference illustrates further connections to Kirman’s ant colony 
model [31] (which has potential relevance to economic markets) and even the kinetic 
Ising model [32] . 

Our focus here is on continuous-time models with dependence on the current over 
the whole history. In the single-particle case we note that even when the particle remains 
stationary, the time-averaged current changes due to the continuous increase of time t 
in the denominator of J jt. The dynamics of the particle can be thought of as a type of 
continuous-time random walk (CTRW) or “semi-Markov” process with a complicated 
non-exponential distribution of waiting times which, in general, also depends on the 
time of the last jump (so that successive waiting times are not identically distributed). 
This correspondence is particularly clear in the case of a random walker moving only 
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in one direction (see Appendix A) and provides a possible route to genuine continuous¬ 


time numerical simulations rather than the brute-force approach of using a discrete-time 
update rule with very short time steps. However, the situation is more complicated for 
many-particle systems, or those with a dependence on multiple currents, and it may be 
practically difficult to obtain explicit forms for the relevant waiting time distributions. 
There is a vast body of work on CTRWs with identically distributed, typically power law, 
waiting times (see [33l|3l] for just a couple of recent examples, discussing different scaling 
regimes and the effects of bias) as well as on more general time-homogeneous semi- 
Markov processes [35] and applications [36]. Helpful explanations of the connections 
between different commonly-employed formulations can be found in [37] and [38] . 

Feedback based on the time-averaged current clearly has the potential to introduce 
long-range temporal correlations and one might ask how these modify the current large 
deviation principle ([1]), if indeed such a relationship still exists. The chief result of [T 6 ] 
was that if, for some 7 , the limit 

1 /■* 

Iw{q)iq + rq') dr ( 2 ) 


I(j) = lim min 

t—>-oo g(r) t'^ 


to 


exists (and is not everywhere zero), then it is the rate function for a modihed large 
deviation principle with speed F. In other words, we now have 

'J 
t 


Prob 


= j ~ e 


-lUF 


(3) 


where 7 is not necessarily equal to unity. Note that, in (|2|), is the Markovian 
rate function evaluated with transition rates w{q), and g(r) is a trajectory in the space 
of time-averaged currents with hxed initial condition (g(fo) = jo) and hnal condition 

{q(t) = j)- 

This result can be derived heuristically by what has been dubbed a “temporal 
additivity principle” (a time-based analogue of the spatial additivity principle of 
Bodineau and Derrida [39]) in which one notes that the time-averaged current changes 
very slowly for large times so can be approximated as constant over time slices long 
compared with the dynamics. Carefully taking the limit t ^ 00 such that both 
the length and the number of the time slices becomes inhnite, this quasistatic (or 
adiabatic) argument gives an integral form for the probability of seeing a given path in 
current space. Furthermore, in the long-time limit a particular current fluctuation is 
overwhelmingly likely to be realised by the optimal (or typical) path which is found by 
minimizing over all g(r) consistent with the required initial and hnal current conditions. 
For further details of this analysis we refer the interested reader to [I6]I3 A proof also 
appears possible at a more rigorous mathematical level by employing older sample path 
large deviation results of Mogul’skii [42] . 

A natural assumption is that the optimal path minimizing the integral in ([2]) is 
arranged so that g(r) is, in some sense, as close as possible to the temporally local mean 


f Note that the technical assumptions involved may break down in models, such as the zero-range 
process, with infinite state space and dynamical phase transitions [40l|41]; particular care should be 
taken in such cases. 
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current j-uj{q), i-e., the expected current for fixed rates w{q). Expanding about 
substituting in ([ 2 ]) we then have 

I 2 


Hj) 


lim min — 

q{T) F 


[q + rq' - ' 

2D., 


dr 


and 


(4) 


where -Diu(q) is the diffusion constant corresponding to rates w{q). This form clearly 
reveals the similarity with the spatial additivity result of |39] but it is worth emphasizing 
that, in the present context, it is only an approximation. For general final current 
j it is impossible to find a minimizing path g(r) which asymptotically converges to 
Indeed this is already obvious for fluctuations far from the mean in the standard 
Markovian case. 

Whilst ([2]) and ([3]) may seem to be a powerful general result, their direct application 
is somewhat limited in practice since, even for those models in which the corresponding 
Markovian rate function is known, the Euler-Lagrange equations involved in the 
minimization are typically too complicated to be solved analytically whether or not 
the form (jlj) is used. The known exceptions [I 6 ] include various types of history- 
dependent random walk, including those where the current for left and right jumps 
is counted separately. One particularly simple case discussed in Appendix A is the 
unidirectional model with rate v{j) = aj which already demonstrates the existence of 
a large deviation principle with 7 smaller than unity for a range of “strong” memory 
dependence (1/2 < a < 1). Physically, this corresponds to a transition to superdiffusive 
behaviour where the fluctuations of integrated current (equivalently, the position of the 
random walker) scale faster than linearly with time. Such a transition was already seen 
in the elephant random walk and related models [711221123]. 

In the next section, we show how these features emerge from a more general 
approximate analysis which involves an expansion about the hxed points of the dynamics 
and can easily be applied to complicated many-particle systems. 


3. Fixed point analysis 


Lightening the notation by dehning /(g) := jw{q), it is intuitively clear that a fixed point 
of the current must obey 

q = f{q)- (5) 


In other words, the expected time-averaged current flowing in the next infinitesimal 
time interval must be the same as that observed in the past. We denote a fixed point 
value satisfying (|5]) by j* and now turn to examine its stability which, as illustrated in 
figure [H is determined by the slope 
df 


A* : = 


dq 


( 6 ) 


9=J 


Specifically, if A* < 1 (left panel of figure [T]) then fluctuations above the fixed point 
yield on average an instantaneous current f{q) which is smaller than the historically- 
averaged current q and thus there is a reduction back towards the fixed point. Similarly, 
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Figure 1. Sketch of current fixed points given by the intersection of the function 
/(?) ■= jw{q) with the diagonal q: stable (left) and unstable (right) cases. 


a fluctuation below the fixed point has /(g) > g so on average the current increases back 
towards the hxed point. Hence, a hxed point with < 1 is stable and, by the reverse 
argument, one with > 1 is unstable (right panel of figure [I]) This heuristic picture, 
which is essentially the continuous-time version of a “cobweb” stability analysis for a 
discrete mapping, can be made more precise by considering the differential equation for 
the time-dependence of the expected current. This latter conhrms that decay towards, 
or growth away from, a fixed point is typically power law in nature which is physically 
related to the fact that the time-averaged current changes more slowly as time increases. 

It is relatively easy to construct models with multiple stable hxed points whose 
selection is inhuenced by the early-time behaviour, cf., e.g., [2Hlll3l[26] for the discrete¬ 
time case. This is especially true in the case of non-monotonic current dependence or 
multiple currentsjJJ] However, in this paper, we specialize to systems in which there is a 
unique stationary state corresponding to a stable hxed point of the dynamics with some 
current j* and slope A* less than unity. In this case we can expand g(r) about j* in the 
numerator and denominator of (0]) and keep terms to leading order to obtain 


Hj) 


lim min — 

t^oo q[r) 


to 


[{l-A*){q- 3 *) + Tq'\ 
2D* 


dr 


(7) 


where D* := = (/"(^.^ (j)\j=j*) ^ is assumed non-zero. Although we are now 

guaranteed to get a Gaussian form for J(j) this approach should correctly capture the 
scaling behaviour of small huctuations and, in particular, the dependence on A*. 


§ This argument implicitly assumes that the system decays to stationarity on a timescale which is 
short compared with the rate of change of the time-averaged current, so that the instantaneous current 
is well described by /(j). This is equivalent to the quasistatic assumption of the temporal additivity 
principle and, at least for finite state space, should always be true for long enough times. 

II For example, a bidirectional continuous-time random walk in which the hopping rates right and left 
depend separately on the time-averaged number of jumps right and left as vrIjrDl) = ujR/{jR + ju) 
and vlUrAl) = ajL/ijR+jr), respectively, has fixed points for and ji satisfying -a < Jr-Jl < a 
(with Jr + Jl = u). It can readily be checked, via exact minimization, that the rate function for the 
net current j = jr — jl is zero for the corresponding range of values. 
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The minimization in ([7]) is straightforwardly carried out; the corresponding Euler- 
Lagrange equations are linear and yield an optimal current path of the form 

q{^) = j* + + K2T^*~^. (8) 

Here the integration constants are determined by the boundary conditions (q'(to) = JO) 
q{t) = j) as 

{jo - - {j - 


Ki = 

K 2 = 

We now substitute 
rt 


t. 


1-2A* 

0 




-2A* 


{jo - J*)tf - {j - J*)t^* 

^2A*-1 _t2A^-l 

JHD into the integrand of (I 
(1 - 2A* 


( 9 ) 


( 10 ) 


and carry out the integration to find 


'to 


Iw{q){q + rq') dr = 


2D* 


-{K,)\t 


2(A-2A* 


t 


l-2A*i 
0 1 


( 11 ) 


Finally, inserting the form ([9]) for Ki reveals that the right-hand side of (ITT]) scales 
asymptotically linearly with t (so we need 7 = 1 for a non-zero limit) for A* < 1/2, and 
(so 7 = 2 — 2A*) for > 1/2. To be precise, we end up with a modihed large 


2-2A* 


as t 


deviation principle of the form 


Prob 


Jt 


= 3 


exp 


exp 


{l-2A*){j-fY^-' 


{2A* 


2D* 

■1)(J 


for A* < 


T \ 2 


2D* 


.2A*-1.2-2A* 
-Iq I 


( 12 ) 


for A* 


> i 


Physically, for A* < 1/2, there is diffusive behaviour with a modihed diffusion coefficient 
D*/{1—2A*). We see clearly here that A* quantihes the effective strength of the feedback 
- for A* negative, huctuations are suppressed whilst, for A* positive, they are enhanced. 
For A* > 1/2, there is super diffusive behaviour which retains an ageing-type dependence 
on the initial time to. At = 1/2 one expects logarithmic corrections whose analysis 
is beyond the scope of the current paper. 

As mentioned earlier, this transition is consistent with that already observed in the 
single-particle example of Appendix A and other random walk models [TJ |22l [23l [26] . 
In the next section we will illustrate the power of the general approach by appeal to a 
specihc many-particle system. 


4. Exclusion process with current-dependent memory 

4 . 1 . Model 

The totally asymmetric simple exclusion process (TASFP) was hrst introduced in 1968 
to describe protein synthesis [H] and, since then, has enjoyed widespread success both 
as a base model for various transport processes ^51 H6] and as a vehicle for advancing 
theoretical understanding of non-equilibrium systems see, e.g., nzmiiiiHi and references 
therein. We here start from the standard continuous-time version of this model on a one¬ 
dimensional lattice with open boundaries and modify it to include a current-dependent 
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Figure 2. Schematic of one-dimensional TASEP with input rate depending on time- 
averaged current j over the whole past history. 


input rate. Obviously, many other forms of current dependence could be envisaged but 
this is a natural choice as a form of feedback - the reader is invited to imagine controlling 
the arrival of cars onto a stretch of road. 

To be more concrete, our model is dehned in the following manner (see also hgure|2]). 
Each of the L lattice sites has only two possible conhgurations: occupied (particle) or 
vacant (hole). A particle on site I hops after an exponentially distributed waiting time 
with mean 1/p to site / -|- 1 if, and only if, that site is vacant. Without loss of generality, 
we set the rate p = 1 in the following. Particles are removed at the right-hand boundary 
(site L) with rate 13 and injected subject to the exclusion rule at the left-hand boundary 
(site 1) with a rate a{j) which, crucially, is a function of the time-averaged input current 
over the whole previous history. In fact, it is obvious from the continuity equation that, 
for a hnite chain in the long-time limit, the time-averaged current must be the same 
between any pair of nearest-neighbour sites. 

For illustrative purposes (and in analogy with the single-particle analysis 
of Appendix AD we mainly consider a linear current dependence of the form 

a(j) = ao + aj, (13) 


where 0 < ao < 1 and a > 0, before later touching on some other choices. An 
input rate of apparently similar form to (IT^ was independently proposed by Sharma 
and Chowdhury [49] to model the recycling of ribosomes in protein synthesis and 
implemented for the more general case of the /-TASEP with extended objects. We 
remark here that in |39| one has the restriction a < 1 (as behts the biological context) 
and also, signihcantly, j is the instantaneous mean (output) current rather than the 
average over the whole previous history. The relevance of these distinctions should 
become apparent in the discussion of phase diagrams and fluctuations below. 


4-2. Mean current 

It is well known (see, e.g., Id) that, in the thermodynamic limit, the standard 
Markovian TASEP has the following three regimes. 

• For a < 1/2, (3 > a there is a low-density (LD) phase in which the mean current is 
controlled by the input rate and given by a{l — a). 

• For a > (3, (3 < 1/2 there is a corresponding high density (HD) phase in which the 
mean current is controlled by the output rate and given by /9(1 — (3). 
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• For a>l/2, /3>l/2 the system is in the maximal current (MC) phase where the 
mean current is limited by the bulk hopping rate and given simply by 1/4. 

We now seek to determine the effect of the current-dependent memory on the phase 
boundaries and the mean current in each phase. 

Following the approach of the previous section, we argue that the mean current in 
the long-time limit is given by the fixed points in the three different regimes: 

' «(/*)(! - «(/*)) for «(/*) < > «(j*) [LD] 

= < for a(j*) >/3,/3 < I [HD] (^^ 4 ) 

^ for a(j*) > i,/3 > i [MCj. 

Unsurprisingly, since the current dependence is in the input rate, the fixed point is 
unchanged in HD and MC phases. In the LD phase, however, some simple algebra 
yields 

— ( 2 q:oO -|- 1 — n) -|- dccgO H (1 — o)^ 


3 


3 = 


2 a? 


(16) 


with the other solution to the quadratic corresponding to an unphysical negative current. 
We can readily show that at the value of j* given by flThll 


d3 


— 1 — -y/4 o^q(Z “f" (1 — Cl) '^ 


(16) 


J=J 


For 0 < ao < 1/2 —a/4 we have 0 < < 1 so this is a stable fixed point with “positive 

feedback. Here the upper bound 
1 a 


CDO — 


4 


(17) 


corresponds to the LD-MC phase transition (determined by a{j) = 1/2). Furthermore 
the LD-HD transition line {P = C({j*)) becomes curved rather than straight and is given 
by 

■(1 — a) -|- -^/daoa -l- (1 — a)^ 


/9 = 


2 a 


(18) 


As might be intuitively expected, the general effect of the positive feedback resulting 
from the aj term is to increase the size of the maximal current phase. However, 
as exemplified by the representative cases in hgure [3l we predict the following three 
qualitatively different forms of phase diagram depending on the value of a. 

• For 0 < a < 1, the phase diagram reproduces that given in [49] - the distinction 
between dependence on historically-averaged and instantaneous current is irrelevant 
for calculation of the fixed point although not for the fluctuations (next subsection). 
The LD-HD transition line always passes through the origin and the phase diagram 
reduces to the Markovian case when a = 0. 


• For 1 < a < 2 there is a qualitative difference in that the LD-HD phase transition 
line intersects the /3 axis at /3 > 0. The feedback is strong enough to ensure a 
non-zero mean-current in the LD phase even for —)■ 0; at cdo = 0 there is an 

unstable fixed point at zero and a stable fixed point at j* = (a — l)/a^. 
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MC 


HD 


Do 

(b) a = 1.6 





HD 

Do 

(c) a = 2.4 


Figure 3. Phase diagrams for current-dependent TASEP with a{j) = ao -I- aj and 
different values of a. Note that the picture in (c) is unchanged for all a > 2. 


• For a > 2 there is no LD phase. In other words, Dq never controls the current - for 
/3 < 1/2 it is determined by the output rate and for > 1/2 by the bulk hopping 
rate. 

The hxed points in the different regimes of these phase diagrams are conhrmed by 
Monte Carlo simulations. For example, in hgure 0] we show a three-dimensional plot 
of the hnal current for a single long trajectory (as a function of boundary rates Dq 
and P) in the model with a = 0.8 and L = 1000; the accord with the theoretically 
predicted phase boundaries is self-evident. As a further check, we then plot in hgure [5] 
the mean time-averaged current from 1000 different histories for the cross-section of the 
phase diagram with f3 = 0.6. The quantitative agreement of the mean current with the 
predicted hxed point j* ffT^ is very good and similar conhrmation is found for other 
rate parameters. However, due to the size of the system, one does need to simulate for 
relatively long times until the rate of change of the current is slow compared to the decay 
to stationarity and the quasistatic assumption is reasonable]^ For smaller systems, the 
decay to stationarity is obviously faster but there are hnite-size corrections for the mean 
currents [50] which would require L-dependent expressions on the right-hand side of (fT4|) 
and, in general, numerical solution for the LD hxed point. For diherent initial conditions 
the current may be diherent for short times but should eventually approach the same 
stable hxed point except for in the special case where the system is started exactly at 
an unstable hxed point. This latter is relevant for simulations at dq = 0 in the a > 1 
case where an initial condition of jo = 0 is observed to lead to a zero current for all 
times whereas jo > 0 gives convergence to the stable hxed point j* = (a — l)/a^. 

In concluding this subsection we note that modihed phase diagrams have been 
calculated for many other variants of the TASEP including those with stochastic gating 
(which can be thought of as the introduction of additional “hidden” variables in the 
standard Markovian model) [51] and density feedback control [52]. However, we stress 
here that our approach enables us not only to predict the mean current but also to gain 
information about the huctuations, as we shall see in the next subsection. 

^ Formally the method requires that the long-time t —> oo limit is taken before the thermodynamic 
L ^ oo limit; this is also important for the study of ffuctuations in the next section. 
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Figure 4. Monte Carlo simulation data for final value of time-averaged current F/t 
as a function of rates ao and /3 for a single trajectory of length t = 10® in a system 
of size L = 1000 with a{j) = uq + 0.8j. Initial condition used was to = 1, jo = 0, 
and each site independently occupied by a particle with probability corresponding to 
the bulk density of a Markovian TASEP with the same input and output rates; for 
times T > to a discrete-time random sequential update rule was used with 20 steps per 
unit time up to r = 1000 (allowing for the fact that the time-averaged current ^(t) 
changes relatively fast at the beginning of the trajectory) and 2 steps per unit time 
thereafter. Data sampled at boundary rate increments of 0.02 and interpolated with 
gnuplot. Solid black lines show theoretically predicted phase boundaries. 


4-3. Fluctuations 


According to the analysis of the Markovian TASEP in [53], the diffusion constant in 
the MC phase scales asymptotically as so, in the thermodynamic limit, D* —)■ 0 

and (I7j) is not applicable. However, in both HD and LD phases the diffusion constant 
approaches a hnite limit - here we aim to understand the effect of the memory on the 
fluctuations in the latter case. 

Starting from the Markovian result in [53] we have, for the LD phase. 


D* = aij*){l-air))il-2aij*)) (19) 

where, for our model with a{j) = ao + aj, the hxed point j* is given by ffT5l) . Now, 
as argued in section [S] we expect long-time diffusive behaviour with modihed diffusion 
coefficient D*/{1 — 2 A*) when A* of (IT^ is less than 1/2. Signihcantly, however there 
should be long-time super diffusive behaviour for 1/2 < A* < 1 which is true for 


ao < 


1/4 - (1 - ay 
4a 


— • Oir 


( 20 ) 
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Figure 5. Monte Carlo simulation data for mean {F)lt as a function of oq with 
a{j) = ao + 0.8j and /? = 0.6. Data from 1000 trajectories of length t = 10® in a 
system of size L = 1000. Other simulation details same as those used for figure 01 
Blue dashed line is theoretical prediction of (fl^ for j*. 


Note that etc is positive only for 1/2 < a < 3/2; in that range, we predict a subregime in 
the LD phase for which the fluctuations are super diffusive and the variance of the time 
averaged current Jjt scales in the long-time limit as This asymptotic scaling is 

fairly convincingly supported by a log-log plot of variance against time for selected values 
of cto in the a = 0.8 case (hgureE]). Additionally, hgure [7] shows a naive check on the 
predicted coefficient D*/\1 —2A*\ across a constant-/9 cross-section of the phase diagram 
with the divergence at etc clearly to be seen. The slight theoretical overestimation of 
the t = 10® data for small cto (corresponding to small mean current) may be related to 
the fact that, for totally asymmetric systems such as this, the current distribution must 
be cut off at j = 0 and the Gaussian approximation is thus expected to be less good for 
small j* (and inapplicable for j* = 0). 

For the same model with a = 1.6, preliminary simulations (not shown) support 
the assertion that there is no superdiffusion and, in fact, suggest that the width of the 
time-averaged current distribution decays somewhat faster than the diffusive prediction 
of (|T2|) . at least for intermediate times. Again this may be related to the j = 0 cut-off 
but further investigation for longer times would certainly be desirable. More generally, 
the existence of a superdiffusive subregime in the phase diagram clearly depends on 
the precise form of the current dependence. For example, in another tractable case 
a{j) = ao + a^/j we also see the MC phase extended at the expense of the LD phase 
but predict that fluctuations throughout the LD phase remain diffusive for all values 
of a. In this case too, the mean current and absence of superdiffusion are conhrmed 



Var(J)/t, Var(J)/t2^* Var(J/t) 
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Figure 6. Variance of J/t as a, function of time for selected values of ao with other 
parameters as in figure [5] (a = 0.8, /? = 0.6). Points show simulation data for (top to 
bottom): uq = 0.01, 0.05, 0.08, 0.12, 2. Black solid lines are fits corresponding to power 
laws with negative exponent min(l,2 — 2A*)\ logarithmic corrections are expected at 
the critical point ac = 0.065625. 



Figure 7. Monte Carlo simulation data for variance oi J as a function of ao for 
parameters of figure [5] (a = 0.8, /3 = 0.6). Red + symbols show 
expected to have finite long-time limit in diffusive regime (a > ac = 0.065625); green 
X symbols show , predicted to be finite in superdiffusive regime 

(a < ac). Blue dashed line is theoretical prediction for D*/\1 — 2A*\ from (ITCl) and (fTOl) . 




Fluctuations in interacting particle systems with memory 


14 


by simulation but more work is still needed to definitively determine the applicability 
of ([12]). 


4-4- Exact numerical minimization 


Going beyond the mean and diffusion coefficient, the current large deviations in the 
Markovian TASEP have remarkably been fully characterized recently for all hopping 
rates and systems sizes m ESI [56]. We can now use these results to evaluate (I2|) 
directly and thus to check the consistency of the Gaussian approximation applied in the 
previous subsections. 

In the LD phase the scaled cumulant generating function (the Legendre transform 
of the rate function) approaches an L-independent limit as the system size increases. 
Specifically, we have 

lim - log(e“^'^) = ail — a) ( - - -for — log ( -< A < oo (21) 

t^oot ^ ’\\-aFae-^ a ^ ’ 


which straightforwardly corresponds to 


U3) 


(2a 


i) + yr^ 

2 


f j log 


(l-a)(l 


- 2j - yT^) 
2 aj 


for 0 < j < 1/4. (22) 


The form (|2TD was obtained by Bethe ansatz in [5l| and via a general parametric 
representation in [55] . It can also be derived within the framework of macroscopic 
fluctuation theory EH. One can readily check that the rate function fl22|) has a zero 
at the mean current ja = a{l — a) with (inverse) second derivative at that point 
corresponding to the diffusion coefficient a(l — a)(l — 2 a). At j = 1/4 there is 
a dynamical phase transition to a regime in which Ia{j) retains a dependence on 
L. Nevertheless, at least away from this transition we claim that the rate function 
of the non-Markovian current-dependent process should be given by minimizing the 
integral in (|2|) with an integrand Ia{q) (<? + rq') which is simply obtained from (12^ via 
the replacement of a with a{q). In practice, this integral is much too complicated to 
approach analytically so we resort to exact numerical calculations using Mathematica. 
Some computational difficulties are encountered here, apparently related to stiffness of 
the differential equations (as well as perhaps the finite range of applicability for Ia{j) 
and the impossibility of negative currents). However, notwithstanding this, the approach 
enables us to push the bounds of investigation beyond the Gaussian regime discussed 
above. 

Returning to our favourite example with a{j) = ao + aj, we focus now on small 
values of a because they lead to more stable numerics and yet clearly illustrate the effect 
of even weak memory dependence on the current large deviations. Figure [H] shows the 
finite-time quantity 


J(j, t) = min - / 4(g) {q + rq') dr 

t Jto 

evaluated at t = 1000 for fixed ao and both zero and non-zero values of a. 


(23) 
In the 
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Figure 8. Mathematica results for LD-phase I{j, 1000) in case a{j) = oq + uj with 
oo = 0.2, a = 0 (black + symbols) and a = 0.1 (red x symbols); initial condition 
to = 1, jo = 0. Black dotted line is exact expression (12^ for Markovian rate function 
/^□(j); red dashed line is Gaussian expansion of non-Markovian /(j) for a = 0.1. Inset 
shows close-up around mean current. 


a = 0 case we find excellent agreement with the Markovian rate fnnction fE2]) and we 
anticipate that onr nnmerical method converges fast towards the long-time limit also in 
the non-Markovian case. For a > 0 the mean current (zero of the rate function) is shifted 
to a higher value and the width of the distribution increased. The approximation from 
the hxed point analysis of the preceding subsections matches very well the behaviour for 
small fluctuations but, as expected, is inaccurate for larger fluctuations. In particular, 
by construction, the Gaussian fails to capture the asymmetry of the rate function about 
the mean - we need the full minimization to see that the probability of large fluctuations 
below the mean is hardly affected by the memory whereas large fluctuations above the 
mean become much more likely than in the Markovian case (presumably because, in 
this model, the feedback can increase but never decrease the hopping rate). 

As a second example, we take a current dependence which illustrates the possibility 
of negative, as well as positive, feedback. Specihcally we set 

a{j) = (24) 

where jao,p is the mean current of a Markovian TASEP with boundary rates ao and (3. 
Note that a{j) thus has a different expression in each of the three regimes of the (ao, (3) 
phase diagram. For any choice of these boundary rates it is easy to see that j* = jao,p 
is a hxed point for all k and, in fact, using the now-established method we hud that 
for K < 8 this hxed point is always stable. In other words, for k < 8 the mean current 
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j 

Figure 9. Mathematica results for LD-phase I{j, 1000) in case a[j) = 
with tto = 0.15, K = 0.5 (red x symbols) and k = —0.5 (green + symbols); initial 
condition to = 1, jo = 0. Black dotted line is exact expression (E21) for IctoU): coloured 
dashed lines are Gaussian expansions of I{j) for k = ±0.5. Inset shows close-up around 
mean current. 

and phase diagram are identical to the nnderlying Markovian model bnt, of course, the 
fluctuations are different. In the LD phase, we have 

A* = fi;ao(l - 2ao) (25) 

and so, for k > 4, there is a super diffusive subregime centred around Oq = 1/4. On the 
other hand, for k < 4 we predict diffusive fluctuations throughout the LD phase with 
modihed effective diffusion coefficient 

D* _ iao(l — Qo)(l — 2 qo) 

1-2A*~ 1 - 2ao(l - 2ao)K ■ 

In accordance with intuition, negative values of n act to suppress fluctuations and reduce 
the width of the distribution about the mean current while positive values promote 
fluctuations and increase the width of the distribution. This is conhrmed by the results 
shown in hgure [9] which again demonstrate that the Gaussian approximation agrees 
closely with the full numerical minimization for small fluctuations but not for large ones 
(especially below the mean). 

5. Discussion 

In this paper we have investigated some aspects of a class of interacting particle systems 
where the rates depend on the time-averaged current J/t. This memory dependence 
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is effectively a form of feedback which can act to suppress or enhance fluctuations. 
In particular, we here considered the application of a recently proposed “temporal 
additivity principle” [16] for obtaining the large deviation rate function for current 
fluctuations in such non-Markovian models via a minimization involving the Markovian 
rate function. Using a heuristic analysis based on hxed points of the dynamics we 
detailed how a Gaussian approximation for the behaviour of small fluctuations emerges 
from the full minimization problem and were thus able to highlight the conditions for 
long-time superdiffusive behaviour. Whilst this approach fails in general for large 
fluctuations, it nevertheless provides a means to gain some information about the 
effects of memory even when the full Markovian rate function is unknown or analytical 
minimization impossible. This claim was corroborated by checking the predictions with 
simulation data for the current mean and variance in a paradigmatic exclusion process 
model, as well as comparing corresponding exact numerical minimization results. In 
order to explore further the underlying assumptions for the temporal additivity principle 
and its approximation, it would be interesting both to put the central arguments of this 
paper on a more rigorous mathematical footing and to develop computational methods 
(perhaps along the lines of the “cloning” algorithm [581159] for Markovian models) to 
efficiently access the full rate function in simulations. 

Although our chief example here was the totally asymmetric simple exclusion 
process one can also apply similar considerations to models with partially asymmetric 
dynamics. Indeed, both the range of applicability of the Gaussian approximation and 
the stability of numerics would potentially be improved without the cut-off at j = 0. 
In the context of jumps in both forward and backward directions, a topical question is 
whether one hnds a so-called fluctuation relation [601 EH E2] governing the probabilities 
of positive and negative currents. Within the Gaussian approximation developed above 
(for A* < 1), we hnd 


Prob(Jt/t = -j) ^ 

exp 

r 2 1 - 2A*)j* 1 

X Jt 

D* 

Prob( J)/t = j) 


D* ° 

exp 


for A* < I 


for A* > i. 


(27) 


This suggests the standard Gallavotti-Gohen-type fluctuation symmetry for A* < 1/2 
and a modihed form for A* > 1/2. The latter is reminiscent of a similar hnding for 
anomalous dynamics in a different setting [63] but a word of caution is necessary here. 
Any Gaussian distribution for the current will necessarily have a ratio between positive 
and negative currents whose exponent is linear in j. This does not guarantee that the 
same symmetry holds in the non-Gaussian tails of the distribution which are neglected 
by this approximation. For time-homogeneous GTRWs or semi-Markov processes (with 
hnite state space and hnite mean waiting time), earlier work [Ml [65] asserts that a 
sufficient condition for the standard symmetry in the full current distribution (arising 
from a time-reversal relation at the level of microscopic trajectories) is the “direction¬ 
time independence” property [381 [35]. In our framework, an exactly solvable model 
with the analogous condition that the ratio of jumps left and right is a constant 
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(see also Appendix A ), was indeed found in (TB] to obey the symmetry. It would 
be interesting to see if there are similar necessary/suflicient conditions for a modihed 
symmetry relation in the case of super diffusive fluctuations. However, in general it is 
not clear that any such relation exists, much less that it has a meaningful interpretation 
in terms of entropy (cf. the discussion in [6l]). 

Other scenarios worthy of closer attention include fluctuations in models with 
multiple stable hxed points and fluctuations beyond dynamical phase transitions. In the 
latter case, one anticipates the possibility of non-convex rate functions corresponding 
to non-differentiable points in the scaled cumulant generating function^ Usually 
for Markovian models, a Maxwell-type construction gives phase separation in time 
and a linear section in the rate function but the introduction of long-range temporal 
correlations means the phase boundary may acquire a hnite probabilistic cost even in the 
long-time limit. This is completely analogous to the manner in which long-range spatial 
correlations can give rise to non-concave entropies in equilibrium [67]. Finally we remark 
that, as already mentioned in [16] the additivity formalism should also be applicable 
to intrinsically non-Markovian models, such as the Alzheimer random walk [68l El |69] , 
but with a non-local minimization problem involving delay differential equations. The 
full analytical treatment of such problems appears an even more formidable task but a 
stability analysis of the dynamics could provide some hope. 

Understanding fluctuations in systems with memory is clearly important from 
both foundational and practical viewpoints but there is much work still to be done in 
establishing connections between different approaches (especially from complementary 
mathematics and physics traditions) as well as in forging new ground. In this context 
we expect that the workhorses of statistical mechanics such as random walk models and 
exclusion processes will continue to play an important role. 
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Appendix A. Details of a single-particle model 

Whilst the bulk of this article is concerned with many-particle systems, we here illustrate 
the formalism by presenting some details of the calculations for a single random walker 

For a mathematical demonstration of a non-convex rate function appearing in another type of non- 
Markovian model, see [66] . 
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on an infinite one-dimensional lattice. We focus in particular on unidirectional dynamics 
where the particle hops in continuous time to the right only with a rate v. In the 
memoryless case the rate function for the number of jumps made by the particle is just 

Iv{j) = V - j + j\n- for j > 0 (A.l) 

V 

which is easily obtained from the limiting behaviour of a Poisson process. Notice that 
this is a convex function with a zero at the mean current = v. We now modify the 
picture by considering a rate v{j) with some functional dependence on the time-averaged 
past current j. 

First, in order to understand the correspondence to the CTRW picture we examine 
the waiting time distribution. We let r be the time of the last jump with a time-averaged 
current immediately afterwards of q (i.e., the particle last jumped to position qr) and 
seek to hnd the cumulative distribution function (cdf) Fq^r{s) of the waiting time s until 
the next jump. In fact, it turns out to be more convenient to study the complementary 
cdf Fq^T-{s) which is just the survival function giving the probability that the particle 
has not jumped up to time s. Now, since the particle’s position does not change, the 
rate at which the next jump takes place depends on time as v^qr/t) where t = r + s. 
Hence we trivially have 


FqAs) 

ds 


= —V 


T + S 


FqAs) 


(A.2) 


with formal solution 


-^g,T(s) = exp y ^ ^ s > 0. (A.3) 

For specihc forms of v{j) one can then calculate the waiting time distribution explicitly. 
For example, specialising to the linear form v{j) = aj -|- b (with a > 0 and & > 0) yields 




T 


r -|- s 


aqr 


^—bs 


for s > 0 


(A.4) 


so that the usual Markovian exponential decay is modihed by a power-law prefactor. 
Note that, in contrast to a standard (time-homogeneous) CTRW, the waiting time 
distribution has an explicit dependence on the last jump time r. 

As an aside, we remark that by an analogous procedure one can calculate the waiting 
times for a bidirectional random walker with right and left rates given respectively by 
vAj) and vAi)- la fact, the distribution of the time between jumps (in any direction) 
is just given by the expression flA.3ll with the replacement of the function v{j) by 
'^Ai) + vAi)- A time-inhomogeneous generalization of the so-called “direction-time 
independence” condition [381ES] that each transition rate can be written as the product 
of an individual waiting-time independent probability and a common factor giving the 
decay of the survival probability, would then seem to require that vAi) AAi) is a j- 
independent constant, i.e., that the two rates have the same functional dependence on 
the current (up to a multiplicative constant). 
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Returning to the unidirectional case, the Euler-Lagrange equation for the optimal 
path g(r) minimizing the integral ([2]) is 

2rg' 


v'{q) 1 




= 0 . 


In the special case v{j) 
straightforwardly solved 


q + rq' q + rq' 

= aj this is a linear differential 
. For a < 1 one finds 


Prob ( ^ = j 


rsj g 




J>0, 


(A.5) 

equation which is 

(A.6) 


whereas for a > 1 there is no stationary state and no large deviation principle. These 
hndings are easily understood within the fixed point analysis of section [3], indeed the 
mean current condition v{j) = j yields a single hxed point at j* = 0 which is stable 
for a < 1 and unstable for a > 1. The Gaussian expansion is not applicable here since 
fluctuations below the fixed point (i.e., with j < 0) are physically impossible. However, 
one can argue directly from (I A. 6 1) that a transition from subdiffusion to super diffusion 
occurs when the slope at the hxed point exceeds 1/2. 

For the general linear dependence of v{j) = aj + b, with 0 < a < 1, there is a 
stable hxed point at j* = &/(! — a) and one can carry out the Gaussian expansion of 
section [3] where, for this unidirectional random walker, D* = j*. Here, there is still 
a transition at a = 1/2 but in the “weaker” memory phase with a < 1/2 the b term 
means that dihusive huctuations (rather than subdihusive) are dominant. The same 
approach can be carried out with more complicated current dependence as long as there 
is a single hxed point around which v{j) has a linear dependence on the past current j. 
For example, with v{j) = a^/j + b (and a > 0, 6 > 0) the hxed point is 

2 


with a slope 


J = 


A* = 


a + y/o? 


+ 46 


(A.7) 


1 + + 46/a2 

which is constrained to be less than 1/2 so the huctuations are always dihusive, albeit 
with modihed dihusion coefficient D */(1 — 2A*). A discussion of similar phenomenology 
in a many-particle model can be found in section 0] of the main text. 
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